clear all;close all;clc
gavgdir = '/home/predatt/anatod/ISI/MEG/GAVG';
cd(gavgdir)
load subjname
load dataset
load ChanSel

rootdir = '/home/predatt/anatod/ISI';
datadir = '/home/predatt/anatod/ISI/raw'; 

for sb = 1:length(subjname)

infofile = strcat(subjname{sb},'_datainfo.mat'); 
subdir = fullfile(rootdir,'MEG',subjname{sb}); 
cd(subdir);
load(infofile)
cd(datadir)

cond = 'all';
cond_trig{1} = [33 34 35 36];


        cfg = [];
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl;
        trial_indx = [];
        for j=1:size(cfg.trl,1)
            if sum(cfg.trl(j,4)==cond_trig{1})
                trial_indx(end+1)=j;
            end
        end
        
   
        cfg = [];
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl(trial_indx,:);
    
   % preprocess 
        cfg = [];   
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl;
        cfg.continuous = 'yes';
        cfg.channel = {'MEG'};
        tmpdata = ft_preprocessing(cfg);        
        
        if sb == 17
            grad = tmpdata.grad;
        end

   % find neighbours
        cfg = [];
        cfg.method = 'distance';
        neighbours = ft_prepare_neighbours(cfg,tmpdata);

	% redo ICA		
        cfg = [];
        cfg.unmixing     = datainfo.comp.topo;
        cfg.topolabel    = datainfo.comp.topolabel;
        comp_data        = ft_componentanalysis(cfg, tmpdata); 
  
    % remove unwanted components
        cfg = [];
        cfg.component = datainfo.comp.artifact; 
        data = ft_rejectcomponent(cfg, comp_data);
        clear comp_data
        
    % fix the gradiometer definitions    
        data.grad = tmpdata.grad;
        clear tmpdata
  
   % preprocess 
        cfg = []; 
        cfg.channel = {'MEG'};
        cfg.demean = 'yes';                      
        cfg.blcwindow = [-inf +inf];          
        cfg.detrend ='yes';                 
        cfg.dftfilter = 'yes';              % notch filter to filter out 50Hz
        cfg.lpfilter = 'no';                % low-pass filter
        cfg.continuous = 'yes'; 

        data = ft_preprocessing(cfg,data);
        
                
   if sb == 17
        
        cfg = [];
        cfg.method = 'template';
        cfg.layout = 'CTF275.lay';
        cfg.neighbours = ft_prepare_neighbours(cfg);
        cfg.missingchannel = {'MRF32'};
        cfg.method = 'nearest';
        [data] = ft_channelrepair(cfg, data);
        
   data.grad = grad;
   end

    % Calculate planar gradient of each trial
        cfg = [];
        cfg.planar = 'yes';
        cfg.planarmethod = 'sincos';
        cfg.trials = 'all';
        cfg.neighbours = neighbours;
        cfg.method = 'template';
        data_planar = ft_megplanar(cfg,data);
        clear data;
 
%% calculate power
         
    cfg.latency = [-0.5 -0.1675-1/1200];
    bsldata_theta = ft_selectdata(cfg,data_planar);
    cfg.latency = [-0.5 -0.2895];
    bsldata_beta = ft_selectdata(cfg,data_planar);

    cfg.latency = [0.5 0.833-1/1200];
    tfrdata_theta_post = ft_selectdata(cfg,data_planar);
    cfg.latency = [0.18 0.3905];
    tfrdata_beta = ft_selectdata(cfg,data_planar);
    cfg.latency = [0.17 0.5025];
    tfrdata_theta_pre = ft_selectdata(cfg,data_planar);
      
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Do TFR calculations using Hanning taper for low frequencies
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        cfg= [];
        cfg.output      =   'pow';
        cfg.channel     =    'MEG';
        cfg.method      =   'mtmfft';
        cfg.keeptrials  =   'yes';
        cfg.keeptapers  =   'no';
        cfg.taper       =   'hanning';
        cfg.tapsmofrq   = 9; %     
        cfg.foi         = 23; %   

        data_tfr_beta        = ft_freqanalysis(cfg,tfrdata_beta);
        bsl_tfr_beta         = ft_freqanalysis(cfg,bsldata_beta); % baseline

        cfg.combinemethod = 'svd';
        fft_beta             = ft_combineplanar([],data_tfr_beta); 
        bsl_beta             = ft_combineplanar([],bsl_tfr_beta); 

        cfg.tapsmofrq        = 6; %     
        cfg.foi              = 1; %   

        data_tfr_theta_pre        = ft_freqanalysis(cfg,tfrdata_theta_pre);
        fft_theta_pre             = ft_combineplanar([],data_tfr_theta_pre); 
        data_tfr_theta_post        = ft_freqanalysis(cfg,tfrdata_theta_post);
        fft_theta_post             = ft_combineplanar([],data_tfr_theta_post); 
        bsl_tfr_theta             = ft_freqanalysis(cfg,bsldata_theta);
        bsl_theta                 = ft_combineplanar([],bsl_tfr_theta); 

        
     theta_bsl = mean(bsl_theta.powspctrm(:,channels_tfr_beta,:),2);
     theta_power_prestim = mean(fft_theta_pre.powspctrm(:,channels_tfr_beta,:),2);
     theta_power_poststim = mean(fft_theta_post.powspctrm(:,channels_tfr_beta,:),2);

     beta_bsl = mean(bsl_beta.powspctrm(:,channels_tfr_beta,:),2);
     beta_power = mean(fft_beta.powspctrm(:,channels_tfr_beta,:),2);
        
     theta_bsl = 10*log10(theta_bsl);
     theta_power_prestim = 10*log10(theta_power_prestim);
     theta_power_poststim = 10*log10(theta_power_poststim);
     beta_bsl = 10*log10(beta_bsl);
     beta_power = 10*log10(beta_power);
       
     beta = beta_power - beta_bsl;
     prestim_theta = theta_power_prestim - theta_bsl; 
     poststim_theta = theta_power_poststim - theta_bsl; 
 
[rho pval]  = partialcorr(poststim_theta,beta,prestim_theta);

partial_corr_and_p(sb,:) = [rho; pval];

disp(sb)
disp(partial_corr_and_p)

clear poststim_theta beta prestim_theta rho pval

end

null_hypothesis = zeros(24,1);
     
disp(median(partial_corr_and_p))

[h p ci stats] = ttest(partial_corr_and_p(:,1),null_hypothesis)
